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Abstract 

In this paper we design an iterative domain decomposition method for free 
boundary problems with nonlinear flux jump condition. Our approach is related 
to damped Newton's methods. The proposed scheme requires, in each iteration, 
the approximation of the flux on (both sides of) the free interface. We present 
a Finite Element implementation of our method. The numerical implementa- 
tion uses harmonically deformed triangulations to inexpensively generate finite 
element meshes in subdomains. We apply our method to a simplified model for 
jet flows in pipes and to a simple magnetohydrodynamics model. Finally, we 
present numerical examples studying the convergence of our scheme. 

1 Introduction 

In this paper, we propose a numerical iterative method for approximating the 
solutions of free boundary problems in two dimensions. Our iterative method 
for free boundary problems is based on Domain Decomposition and damped 
Newton's method ideas. In general terms, free boundary problems seek to de- 
termine unknown function u with some prescribed conditions on a unknown 
interior interface, exterior boundary or (sub)domain. In many applications, it 
is prescribed the value of u on the free interface and it is required that u sat- 
isfy a condition involving (both sides) derivatives of u on the interface. We 
mention jump conditions of Stefan, Bernoulli and Gibbs-Thomson type, among 
others. There is a considerable literature of iterative methods for these type 
of free boundary problems; see for instance [H [7j HI [13l HH [20j [22] and refer- 
ences therein. In particular, numerical finite elements methods have been pro- 
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posed to solve Stefan-like free boundary problems (including time dependent 
problems) and some other similar phase transition problems; see for instance 
[21 EH HS1 HH HH] . These methods use a variational formulation of their original 
problem. Level set approach for Stefan problems were also proposed in [5] and 
references therein. 

The free boundary conditions that we deal with, up to our knowledge, have 
not been extensively studied from the numerical point of view. We are par- 
ticularly interested in free boundary problems were the unknown function sat- 
isfy nonlinear jump constraints across the free interface. More precisely, given 
fl C R 2 , g : dfl — >• M and A 6 M, we want to find a function u : — > M and a free 
interface T (diving fi in two subdomains fi + = {u > 0} and fi~ = {u < 0}) such 
that u and V satisfy the following subdomain equations and boundary condition, 

-div(a x V«) =0 in Q+, (1) 

-div(a 2 Vw) =0 in O", (2) 

u = g on <9f2 (3) 

and free interface condition 

ai|Vu+| 2 - a 2 |Vu~| 2 = A on V, (4) 

or, similar nonlinear constraint for the jump in the derivative of u across the free 
interface T. Above, u + and u~ denote the value of the solution u on both sides of 
the free interface and the derivatives and the quantities involved are interpreted 
as side limits. In many applications, the interface conditions are imposed in a 
weak sense. These conditions can be also interpreted if we replace the operators 
involved (e.g., trace of the derivatives) by some smooth or regularize version of 
them when necessary. 

We are not aware of a simple inexpensive numerical method to solve prob- 
lem (|T|)- ([4j) . The finite element methods mentioned earlier to handle Stefan, 
Bernoulli and similar free boundary conditions are based on variational formu- 
lations. They do not seem to be easily extended to handle our nonlinear free 
boundary constraint. Also, Bernoulli type free boundary problems when one 
of the phases is a constant function seems to be easier to handle numerically. 
In this case, using the fact that the tangential derivative on the free interface 
is zero and that the flux sign can be a priori determined, the interface condi- 
tion reduces to a linear condition of the form d^u = A where d v is the normal 
derivative on the free interface. 

We have two main applications in mind: 1) the jet flow model studied by Alt, 
Caffarelli, Friedman [2j[TJ and 2) a free boundary problem arising in magnetohy- 
drodynamics studied in [101 112] . These applications are simplified mathematical 
versions of complicated flow models and they focus in the main modeling as- 
pects. Despite of the mathematical simplifications, in either case, the resulting 
model problem above is still complex and finding and understanding solutions 
requires numerical methods. The methods used for this problems should be 
inexpensive and simple. The method presented here is designed having these 
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considerations into account. It can also be easily extended to handle different 
free boundary problems such as the stationary solutions of the Stefan's problem, 
and other similar problems. 

The iterative method proposed in this paper for problem (HJ)- OJ) is based on 
the following simple ideas. Assume the solution u is sufficiently regular, and let 
T denote the free boundary of problem |T]). Since u — on T, Vu + = d v u + r] 
on r, where 77 is the outer normal vector of the region defined by the support 
of u + . Hence, the free boundary condition Q reads 

ai\d v u+\ 2 -a 2 \d v u-\ 2 = A on Y. (5) 

Next, assume we have an approximation T of T dividing £1 in two different re- 
gions fl + and f2~. We also assume that il + and f2~ are connected subdomains, 
and dil + n d£l = S + and <9£1~ n dil = E — . In order to construct an approx- 
imation u of u, we can solve Dirichlet problems (JTJ)- (J3j) in the approximated 
subdomains with homogeneous Dirichlet boundary condition on the approxi- 
mated free interface Y. The solution of these two independent problems give u + 
and u~ . We observe that we do not expect the function u to satisfy condition 
([5]), since Y is only an approximation of Y. Finally, we update the approximation 
of the free boundary by using the quantity a = ai\d v u + (x)\ 2 — a2\d v u~ (x)\ 2 — A 
and a perturbation of Y in its normal direction r/(x). More specifically, we lo- 
cally move r in the direction of f]{x) by a magnitude rcr where r is a positive 
damping parameter. Once the new approximation of Y is obtained we restart 
this procedure. 

The rest of the paper is organized as follows. In Section [2] we describe our 
iterative scheme. Section [3] describes the finite element implementation of our 
method. In Section|4]we present the jet flow model proposed by Alt, Caffarelli, 
Friedman and some numerical solutions for this problem. Numerical experi- 
ments for the magnetohydrodynamics problem studied in |101 112] are presented 
in Section O Section [6] presents some numerical experiments where we study 
convergence properties of our scheme. Finally, we present our conclusions and 
comments in Section [7J 



2 Model problem and iterative method for the 
free interface 

In order to simplify the presentation and fix ideas, we consider the two dimen- 
sional case f2 C R 2 , and the following free boundary problem 



-V • (aiVu+) = 
-V • (a 2 Vu ) = 
ai(d v u + ) 2 - a 2 (d v u~) 2 = 
u = 



in {u > 0} 
in {u < 0} 






A on r = {u 
g on dil. 



0} 



(6) 



We also assume there exist two connected curves S + , S , such that <9f2 = 
£+ U and g\%+ > and g\%- < 0. This model problem, or similar system 
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of equations, appear in different applications. 

We approximate the solution of problem ([6J by constructing a sequence of 
approximations of the free boundary T. Assume we have an approximation of 
the free boundary T, then we solve two independent elliptic problems and use 
the condition ai(d v u + ) 2 — a2(d v u~) = A to updated the approximation of the 
free boundary as follows. 

Assume T n is an approximation of the free boundary dividing the domain 
into two subdomains, f2+ (enclosed by S + ur„) and Q~ (enclosed by £~ ur„). 
We define the n-th approximation of u as follows. In f2+ the function u n solves, 

-V • (aiVttn) =0 in fl+ 



u„ 



g on £+ (7) 



u n = on T r 



In Q n the function u n solves, 



-V • (ot2 Vit n ) = in Q~ 

u n = g on £~ (8) 
u„, = on r„. 

The main idea to define the updated approximation of the free boundary 
r n +i is very simple. First, we define 

a = ai(9 t? u+) 2 - a 2 (d ri u~) 2 - A. (9) 

Here, we use the notation d v u^ as the outward normal derivative of with 
respect to the region f2^. Next, if for instance, a(x) > for some point i £ T,,. 
then we would like to locally update T n such that cr(x) is closer to zero. This 
can be done by decreasing the flux of u n in fl^ and/or increasing the flux of 
u n in 0~ in a neighborhood of that point. We expect to obtain this by locally 
moving the free interface T n in the normal direction outward to fi+ . We define 
the new approximation of the free interface by 

r n +i = {x + T &Vr+ i w i tn x e r„}. (10) 

Here r = r(a) is a small positive parameter, and 7f r + represents the unitary 
normal vector of r„ outward to f2+ . 

Finally, we observe that there are several ways to define Tq dividing the 
domain Q in two parts as desired. For instance, we can take Tq as the zero level 
set of any regular extension of the boundary data g. 

Remark 1 We note that we need only an approximation of a (which requires 
only approximation of the flux). This is important in case u is not regular 
enough to allow the computation of the square of the flux. 

Remark 2 We mention that in J^j it is proved that the solution of problem 
{I]]-^p, in the case a\ = a 2 = 1, is a minimizer of the following functional 
J(v) = J n l v>0 (|Vu| 2 + A 2 ) + l„<o (|Vw| 2 + A 2 ,) dx where A = A 2 - \\, and 
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l t ,>o is the characteristic function of the set [{v > 0} := {x 6 Q,v(x) > 0}], 
(similar for 1 v< q). We have also developed a method for this problem based 
on the minimization of this functional. This was performed by, first, introduc- 
ing a regularized approximation J c of the functional J. Next, we looked for a 
minimum of the functional J e by solving the steepest descent evolution PDE as- 
sociated to this functional. However, the observed numerical results were not 
satisfactory. We also observe that this method requires to solve a nonlinear 
problem for each time step resulting in more computational work compared to 
our iterative method. 

We also observe that our method to solve problem (|6]) can also handle dif- 
ferent problems. For instance, the same ideas apply to the following abstract 
free boundary problem. Let C + and £_ represent two second order elliptic op- 
erators. Assume il C R n , and let <j> and if) : M — > K be two increasing functions. 
Consider the problem of finding u and a free interface T such that 




in 

in ft~ , 

ip(d v u~) = X(x) on the free boundary T 
on d£l. 



(11) 



Here T represents the free boundary separating the two phases, d n Ui represents 
the outward normal derivative with respect to the i-th phase. Up to our knowl- 
edge there is no rigorous studies of such general class of problems. We mention 
that these problems include the stationary solutions of two phase Stefan prob- 
lem (see [TTJ [9l [18j [21] and references therein); and other problems involving 
nonlinear free boundary conditions (see j!7j). 



3 Finite element implementation 

Now we describe the finite element implementation of our iterative method. In 
each iteration we have to approximate the solutions of problems (J7J and © as 
well as a in ©. 

Let n £ N be our iteration parameter, 7^ be a triangulation of ft with 
nodes {x% Af=i and edges {e n ^}^f 1: and let T 1 ^ be an approximation of the 

free boundary T, such that C uf^enj. Here we also assume divides ft 
into two subdomains, (enclosed by E+ur^) and ft~ h (enclosed by S~ur^). 

Set V+' h = {ve V l {T£, n+ h ); v = on dSl+ h }, where P^Tj?, fi+ h ) represents 

the set of continuous piecewise linear functions on T n (the space V$' is defined 
similarly). The n-th approximation of the solution u of ([5]), denoted by u\, 
solves the finite element problems, 

J Q + aiVu^Vz' 1 dx = for all z h G V^ Jl 

u 'n( x ' 1 ) = g( xh ) for a11 xh g s+ > and x ' 1 e Xi ( 12 ) 

u h (x h ) = for all x h G r£, and x h G T£, 
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and 

J n - aaVu^Vz' 1 dx = for all z h G Vq' H 

u^{x h ) = g(x h ) for all x h G and x h e 7^ ( 13 ) 
= for all x h G r£, and x h eT*. 

We define p^" t as an appropriate piecewise linear approximation of the flux 
of h across Fjj; see Appendix^ Analogously, we define p~ h as the discrete 

flux of u^Iq- across T^- We note that, 

n , h 

Pn.h = J2 a tj^ and Pn,h = J2 a nJ^' 
x h er /i x h e pfi 

where the function represents a basis for the space V 1 {Tm^) restricted to 
Define, 

-n= E («i) 2 -K,i) 2 -A 2 )^i. (14) 

The new approximation of the free interface is given by the piecewise linear 
curve, 

r h n+1 = {x + Ta*(x)rf n (x) ; with x G T h n }. (15) 

Here fj r + represents an approximation of the unitary normal vector of Tjj out- 
ward to More specifically, since Tjj is piecewise linear, its normal vector 
fj T h(x) is not well define when a; is a vertices of T„- Different strategies can 
be used to handle this problem, for instance, we can define the normal vector 
as the average of the two adjacent normal vectors of x\ or we can interpolate 
the vertices of r„ by a smooth curve and define the normal vector of T 1 ^ at x 
as the normal vector of the smooth interpolation of T^. In our numerics, we 
implemented the first strategy. 

Next, the triangulation T r [ l +1 is defined such that T^ +1 is the union of edges 
in Tn + i- More precisely, we obtain 7^ l + i from T£ using a harmonic extension 
of the displacement Tcr h ff r h as follows. First, we introduce the vector function 

= (wi,u>%) where each component satisfies 

/„+ ai V^Vz h dx = for all z h G V^Tn, fi+ h ) 

wf(x h ) = for all x h G £+, and x h G T£ (16) 

w}(x h ) = T<7 h (x h )(rj rh (x h )) ■ ej for all x h G T„, and x h G C\T% , 



and 



/„- h a 2 Vuj!}Wz h dx = for all z h G V 1 {T n , Sl~ h ) 

wf(x h ) = for all x h G E — , and x h G (17) 

(x' 1 ) = ro- ft (a;' l )(^h (a;' 1 )) • e} for all x h G T„, and a;' 1 G T£, 
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where e± — (1, 0) and ei = (0, 1). Then we define the nodes of the new triangu- 
lation 

X n+l,j =<j (18) 

The edges and triangles structures of Tj? +1 is inherit directly from ■ 

Finally, we observe that the initial triangulation 7o can be defined as any 
regular triangulation of fl containing vertices on the initial approximation Tq of 

r. 

3.1 Summary of the iterative method 

We now summarize the proposed iteration for a given a tolerance etoi- 

Input: Domain f2 and boundary condition g. 

Output: Free interface approximation, F„, and approximation of 

the solution, u^. 

1. Set up To (and the positive and negative subsets fig" and 

n -). 

2. For n = 1, 2, ... , until convergence, do 

(a) Compute u% by solving (fP2|) and (|T3|) . 

(b) Compute a% in (THl) . 

(c) Compute the triangulation displacement by solving 
(pi) and HIl). 

(d) Set up the new free interface approximation rjj +1 in 
(fT5| and triangulation 7^ +1 in (|Tg|) . 

Here, the convergence criteria is given by ||0^||L°°(r) < e to i- 

3.2 The parameter r 

In order to get some insight on the role of the parameter r we may compare our 
method with a regular Newton's method to solve o{z h ) = 0, where z h represents 
the coordinates of the vertices of the partition belonging to the free boundary. 
Formally, a Newton's method for our problem would consists of the following 
the iteration 

V h z a(z n ' h )(z n+1 ' h - z n - h ) = -a(z n ' h ). 

Here V^cr(z n ' ) represents a formal derivative operator of a with respect of z n ' . 
Assuming it is possible to invert the operator W^a(z n ' h ) we would have 

( z n+i,h _ z n,h^ = _(v^ (7 (z n <' l ))- 1 ( 7(2 n <' 1 ). 
From (TTUj) we conclude that our method satisfies 

z n+l,h _ z n,h = T(T ( z «)^ rn 
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Finally, assuming it is correct to update the z n ' h by moving it toward the 
average of normal directions of adjacent to z n ' h , we expect to obtain a 
damped Newton method by choosing r sufficiently small. 

In our numerics we observed that the parameter r should be chosen suffi- 
ciently small to avoid big variations of the triangulation with respect to in a 
single step. 

4 Applications to jets of two fluids in a pipe 

In this section we apply our method to a simplified version of the jet problem 
for two fluids. We consider a version of the model discussed in pQ. There, it 
is considered the model of two planar flows along an infinite pipe with one free 
interface. Here we use our method to computed approximate solutions in a 
truncated pipe model with some given inflow/outflow data. 

The model for jet flow studied in pQ is the following. Let u denote the 
stream function associated to the irrotational flow of two ideal fluids. The 
regions occupied by each different fluid are represented by the support of u + 
and u~, where u + (u~) denotes the positive (negative) part of u. Let N\ : R — > 
(01,02), with < ci < C2 be a continuous and piecewise C 2 function, satisfying 
lim^oo JVi(tf) = B, !™{Nx{y) - B) 2 dy < oo, and / p °° N[(y) 2 dy < oo. The 
Alt et. all. model assumes the two fluids occupy an infinity semi-strip region 
enclosed by the graph {(Ni(y),y); y > a, a < 0}, and the lines y — a and 
x = — 1. The fluids enter the region at the boundary {y = a}, and the two 
fluids are separated from each other in {y < 0} by a given continuous and 
piecewise C 2 curve N 2 : [a, 0] — > (—1 + S, c 2 ), satisfying dist(N 2 ,N 1 ) > 0. A 
special truncated case of this configuration is shown in Figure [1] (left). The 
problem consists in finding the free boundary separating the two fluids in the 
region y > 0, assuming each flow has constant speed when y — ► oo. More 
specifically, we look for u and A satisfying 

Au = in each fluid 

|Vu + | 2 — |Vm _ | 2 = A on the free boundary separating the two fluids , Q . 
u = Q on {x = —1} and u = — 1 on N\ 
u = g on [-l,iVi(o)] x {y = a} 

where 

f A = 1/(1 + b) 2 - Q 2 /{B - b)' 2 , and the free boundary , . 

1 approximates the point (b, 0) when y — > oo. 

Here the function g £ C 1 is monotone decreasing and 

< g{x) < Q for x < N 2 (a), 
-1 < g(x) < for x > N 2 (a), 
<?(-!) = Q and g(N 1 (a)) = -l. 

Existence and uniqueness of solution for this problem was studied in [2] , where 
it was proved that minimizers of an appropriate functional are weak solutions 
of problem (|T5)) . 
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We construct approximated solutions of the above problem. In particular 
we work with a truncated domain to represent a pipe. 




Figure 1: Simple vertical pipe configuration. 

We now refer to the problem configuration in Figure [TJ Given a positive 
constant Q, and functions N : [i?~,0] — > (—1, 1) and g : [—1, 1] — > K, we want 
to find u : (—1, 1) x (i?~, R + ) — » K and free interface T represented by 

T = {(yJr(y)); with < y < R+} 

where f r : [0,i?+] — > M is such that / r (0) = iV(0) (see Figure [T] right picture). 
The function u and the free interface T satisfy 

Am = in and Au = in fi+ (21) 

where 

0+ = {-1 < x < N{y), R- < y < 0} U {-1< x < fr(y), 0<y<R-} 
fi- = {iV(y) < x < 1, R- < y < 0} U {/ r (y) <a;<l, 0<y < R~}. 

The function u has to satisfy the following known given data, 



u(N(y),y) 


= 0, R 


<y<0; 


(22) 


u{x, Rr) = 




-1 < x < 1; 


(23) 








(24) 


—u(x,R^ 
dy 


) = o, 


-1 < x < 1; 


u {- l ,y) = 


Q, Br 


< y < 


(25) 


u{l,y) = - 


1, R~ 


<y<R + ; 


(26) 



and the following conditions on the free interface 

u = 0, onT (or u(f r (y),y) = 0,0 < y < R + ) (27) 
{d ri u + f - (d ri u-) 2 = A on T (28) 

where A is given by 
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We note that the boundary condition on the top of the domain (see Figure 
[T]) is the homogeneous Neumann boundary condition, hence, the free boundary 
is not fixed at the top. That is, the value fr{R+) is not prescribed. 

For each value of b G (—1,1), we can compute A through (|2TJ)) and use our 
method to find an approximation of u = u b and r = T b (represented by /r = /r) 
that solve ([2~l} - ([2"9")) . Since the free boundary must have the vertical line x — b 
as an asymptote, a feasible approximation of the free interface Y is obtained if 
b = b* where 

(30) 



b* = /* (R+). 

This is compatible with the asymptote condition lim^oo fr(y) = b* . 

Next, we use a bisection algorithm (applied to the function (—1, 1) 3 b t— > 
/p(i? + ) e (—1, 1)) to find the correct value of b* such that ([30| is satisfied. 



In the two examples presented next we run our method described in Subsec- 
tion ED until ||ct!j||loc < tol = icr 6 . 










if 



Figure 2: Free interface that solves (|2T] t - ([29 ]) and {2Q) with N(y) = 0.5|y|/i?~, 
g{x) = {x- 0.5)/1.5 if -1 < x < 0.5, g{x) = 5(x - 0.5)/0.5 if 0.5 < x < 1, 
Q = 5 and b* = 1/3. Initial configuration (Left), final configuration (Middle) 
and a zoom around (0, 0) showing the free interface and the final mesh (Right). 




Figure 3: Function u and free interface Y that solve (|l?Tj) - (l2"!)t and (I3TJ1) with 
N{y) = 0.5(|y|/ir) - 25 , g{x) = (x - 0.5)/1.5 if -1 < x < 0.5 , g(x) = (x - 
0.5)/0.5 if 0.5 < x < 1, Q = 1 and b* = 1. Solution (Left) and zoom around 
(0,0) showing the free interface and the final mesh (Right) 

The first example considers the nozzle represented by N(y) — 0.5|y|/i?~ with 
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Q = 5. The data on the bottom is given by, h(x) = (x — 0.5)/1.5 if — 1 < x < 0.5 
and h(x) = 5(x - 0.5)/0.5 if 0.5 < x < 1. We obtain 6* = 1/3. The initial free 
boundary approximation Tq is the strait line from (0, 0) to (0, R + ) The resulting 
free boundary is displayed in Figure [5J 

In the second example of jet flow problem we consider the nozzle represented 
by N(y) = 0.5(|y|/i?-) ' 25 with Q = 1 and the Dirichlet data on the bottom 
side given by h(x) = (x - 0.5)/1.5 if — 1 < x < 0.5 and h(x) = {x — 0.5)/0.5 
if 0.5 < x < 1. We obtained 6* = 1. The initial free boundary approximation 
To is the strait line from (0,0) to (0,R + ). The resulting free boundary and 
numerical solution (with constant lines -stream lines) are displayed in Figure [3] 



5 Application to a free boundary problem aris- 
ing in magnetohydrodynamics 

In this section we apply our methodology to the model of plasma problem 
studied in [101 112] . Here we are interesting in modeling the plasma confined in 
a Tokamac machine. More specifically, given f2 C R 2 and the positive constants 
7 and A, the plasma problem is to find u, a closed curve T lying in fl and a 
positive constant /3 such that 

—Ait = flu in n~ = int{x £ Q,u(x) < 0}, 
In- u2 = 1 

-Au = in fl + — {x e fl, u(x) > 0}, , . 

u = 7 on dfl, 

u = on r and V = dft~, 

\d v u+\ 2 - \d n u-\ 2 = X on T. 

Here the plasma is enclosed by the curve T, and the complement of this region 
with respect to fl is vacuum. The function u represents a flux function associated 
to the magnetic induction B, satisfying B = (u X2 , —u Xl ,0). 

It is easy to modify our method and apply it to this problem. We follow 
the description in Subsection 13.11 and iterate until ||ct^||l~ < tol = 10 -6 . In 
this problem, the free boundary is a closed curve separating the domain in two 
connected components; as showed in [TU]. The adaptation of our scheme to treat 
this problem is straightforward. We also mention that other formulations of the 
model, having the nonlinear condition on the free boundary, are also possible; 
see [10] and references therein. 

In the first example, we consider the case of f2 being the ball with center (0, 0) 
and radius 1. We choose 7 = 1 and A = 2 2 — l 2 = 3. The initial approximation of 
the free boundary is an ellipse centered at (1/5, 1/5) and with axis 1/3 and 1/2. 
The resulting configuration is depicted in Figure0]and j3 = 13.6727. We observe 
that the final shape of the free boundary approximates a circular region. This 
coincide with the results in [12] where the authors proved that for the domain 
Q being the unit circle, the resulting free boundary is circular and centered 
at (0,0). We note that in this example, the initial approximation of the free 
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interface is far-off from the solution and despite of this fact, our algorithm still 
converges to the solution. 




Figure 4: Free interface that solves ([3T]) with fi = {(x, y) : x 2 + y 2 < 1} ,7 = 1 
and A = 2 2 — l 2 = 3. Initial configuration: an ellipse centered at (1/5, 1/5) and 
with axis 1/3 and 1/2. (left). Final configuration (right). 




Figure 5: Free interface that solves (|3Tj) with depicted domain f2, 7 = 1 and 
A = 5 2 — l 2 = 4. Initial configuration: the circle with center (0,0) and radius 
1/3 (Left). Final configuration (Center). Solution (Right) 

The second example considers the configuration described in Figure [5] The 
domain corresponds to a circle from which it have been cut-off the region 
{y < —2/3} and the intersection with the circle with center (5/3,0) and radius 
1. In this example we use A = 5 2 — l 2 = 4. The initial approximation of the 
free boundary is a ball with center (0,0) and radius 1/3. The resulting free 
boundary is presented in Figure [5] (center) and the solution is plotted in Figure 
(right). The computed value of /3 = 13.7034. 

6 Additional numerical examples 

In this section we present some representative numerical examples. We run our 
method described in Subsection 13.11 until 

\\<t*\\ l °° < toi = 10- 6 . 
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6.1 A known free interface and error decay 



We consider problem (jSJ) with ai = 0,% = 1 and a known exact solution, what 
allows us to measure the accuracy of our method. The domain is fl = [0, 1] x 
[0,1], A = — 1, the boundary data is given by g(x,y) = 2min{x — 0.5,0} + 
max{x — 0.5, 0}. We note that u(x, y) = 2 min{x — 0.5, 0} + maxfa; — 0.5, 0} is a 
solution of the problem. For this exact solution, the free interface is the strait 
vertical line from (0.5, 0) to (0.5, 1). 

We apply our method with the initial approximation of the free boundary 
given by Tq — {(x, y);y — 0.5 + 0.1 sin(27ry)} and the parameter r = 10 -4 . We 
present the initial and final subdomain configuration in Figure [B] In Figure [7] 
we present the L°° norm of log(|cr' l |) in (| 14[) along the number of iterations n. 
We observe a decay of the value \a h \ faster than O(e~°' 004n ). This example also 
show that our stopping criteria is effective in the sense that at the last iteration, 
we see that the L°° norm is already in stagnated plateau for the corresponding 
mesh size. 




Figure 6: Results for the test problem in Section l6Tl Initial configuration (left). 
Final configuration (right). 




Figure 7: Result for the test problem in Section ETTl We plot log(||cr' l || o) over 
the iterations requires to achieve ||<r ||oc < 10 -4 for three different meshes (left). 
Corresponding solution error (in log scale) in | • \h 1 (d) and L°°(D) norms. 
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6.2 An example with heterogeneous coefficients 



This example considers a problem of type © with heterogeneous coefficient in 
each side of the free interface. We consider f2 = {(x,y) : x 2 + y 2 < 1} and the 
coefficient (s) 



The Dirichlet data around the circle is given by g(x, y) = x and we use A = — 1. 
The initial approximation of the free boundary is the strait line Tq = {(0, y), < 
y < 1}. We run our method with r = 10 -5 . We show the resulting free boundary 
in Figure [5] 




Figure 8: Results for the test problem in Section ROl that involves heterogeneous 
coefficients. Initial configuration (left). Final configuration (right) 



7 Conclusions and comments 

We have proposed a simple iterative method to handle free boundary problems 
involving nonlinear flux conditions. It is important to note, that the numeri- 
cal treatment of nonlinear flux conditions on the free interface have not been 
extensively studied in the literature. This is the case despite of the fact that 
the mathematical analysis of simple models with nonlinear flux conditions on 
the free interface have been carried out by Caffarelli and coauthors a couple of 
decades ago. The proposed method is a simple domain decomposition method 
with inexpensive iterations. As a consequence it can be used for the better 
understanding of simplified models of complex flow problems. We present nu- 
merical results showing that, our iterative method is effective and perform well 
in several applications where nonlinear flux jump constrain drive the free inter- 
face behavior. 

We obtained encouraging numerical results with our method but, its math- 
ematical analysis is still needed. In a future work we plan to address mathe- 
matically questions related to the converge of the method. Other interesting 
numerical aspects we want to address are related to the implementation of 
adaptive refinement, the use of inexact local solvers (instead of exact subdo- 
main solvers) , and the design of preconditioners for our scheme. The extension 
to three dimensions can be considered. 




and 




y >o 

y<0. 
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We note that we consider simplified models of complicated flow problems. 
If we want to extend our method for more realistic models we need to consider 
time dependent problems. In this case, it would be important to be able to 
handle topological changes in the evolution of the free boundary. 
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A An approximation of the flux 

Given a free interface approximation T h , we consider the approximation of the 
flux of u h (the solution of problem H2])) on {it. 

Denote by A = [dij] the Neumann finite element matrix defined by 

ai V^V^j dx 

h 

where {&} are the usual hat basis function of the space V 1 (T£, 

We classify the nodes in interior nodes I, boundary nodes S + and interface 
notes r. This classification gives the following block structure of the matrix A, 




A 



£+£ + 



A 
A T 




The solution of (fl"2"|) is given by, 



i$ \ / A n \A^ + g h ) 



We define n by 

/i = Aj r m = AjrAjlA^+g' 1 . 

Let JVr be the number of vertices of on T h . We note that, using basic 
finite element analysis, we see that = (fii) € l Wr with 



(a{\7u h ) ■ dx — I {aiVu h ) ■ rfphcp^ ds. 



Here given i £ {1, ...Np}, £% represents the index of the a node of T h belonging 

to r'\ 
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We use /i to obtain a piecewise linear approximation of the flux Wu h ■ r/ r h . 
Since u h = on T h , for each edge of Ck of T h we have 

Vu h \ ek = d m ur] k 

where rjk represents the normal vector to edge pointing in the outward direc- 
tion of Hence, 

Mi = / {Vrhaiij r h)d V h u(()e i ds. (32) 

We define the piecewise linear approximation of d-q h u as follows. First, 
we observe that A^ G span{(f>e. i \Yh}i<.i<Nr' Next we introduce the matrix Q = 
[qij] e M JVrXJVr with 



9»j = / ( 7 ?rfc £l i J 7r'*)^i^ rfs - 



Finally, based on relation (|3"2"j) we define 

Af - 2 Ir* ( 33 ) 

i 

where a = («j) is the solution of 

= ji. 

In a similar way we define Af, the approximation of of the flux on T h , of the 
solution of lfl3l). 



Remark 3 j4 more regular approximation of the flux can be done in practice. 
For instance, we could obtain a as the solution of the following problem 

(Q + eD)a = f i. 

where D is diffusion of operator on T h and e is a regularization parameter. 



References 

[1] Hans Wilhelm Alt, Luis A. Caffarelli, and Avner Friedman. Jets with two 
fluids, i. one free boundary. Indiana Univ. Math. J., 33(2):213-247, 1984. 

[2] Hans Wilhelm Alt, Luis A. Caffarelli, and Avner Friedman. Abrupt and 
smooth separation of free boundaries in flow problems. Ann. Scuola Norm. 
Sup. Pisa CI. Sci. (4), 12(1):137-172, 1985. 

[3] John W. Barrett and Charles M. Elliott. Fixed mesh finite element approxi- 
mations to a free boundary problem for an elliptic equation with an oblique 
derivative boundary condition. Comput. Math. Appi, ll(4):335-345, 1985. 



16 



[4] F. Bouchon, S. Clain, and R. Touzani. Numerical solution of the free 
boundary Bernoulli problem using a level set formulation. Comput. Methods 
Appl. Mech. Engrg., 194(36-38) :3934-3948, 2005. 

[5] S. Chen, B. Merriman, S. Oshcr, and P. Smereka. A simple level set method 
for solving Stefan problems. J. Comput. Phys., 135(l):8-29, 1997. 

[6] Zhiming Chen, Tsimin Shih, and Xingye Yue. Numerical methods for Ste- 
fan problems with prescribed convection and nonlinear flux. IMA J. Numer. 
Anal, 20(l):81-98, 2000. 

[7] Karsten Eppler and Helmut Harbrecht. Efficient treatment of stationary 
free boundary problems. Appl. Numer. Math., 56(10-11):1326-1339, 2006. 

[8] M. Flucher and M. Rumpf. Bernoulli's free-boundary problem, qualitative 
theory and numerical approximation. J. Reine Angew. Math., 486:165-204, 
1997. 

[9] Avner Friedman. The Stefan problem in several space variables. Trans. 
Amer. Math. Soc, 133:51-87, 1968. 

[10] Avner Friedman and Yong Liu. A free boundary problem arising in mag- 
netohydrodynamic system. Ann. Scuola Norm. Sup. Pisa CI. Sci. (4), 
22(3):375-448, 1995. 

[11] S. L. Kamcnomostskaja. On Stefan's problem. Mat. Sb. (N.S.), 53 (95):489- 
514, 1961. 

[12] Kyung-Kcun Kang, June-Yub Lee, and Jin Keun Seo. Identification of a free 
boundary arising in a magnetohydrodynamics system. Inverse Problems, 
13(5):1301-1309, 1997. 

[13] Kari T. Karkkainen and Timo Tiihonen. Free surfaces: shape sensitivity 
analysis and numerical methods. Internat. J. Numer. Methods Engrg., 
44(8):1079-1098, 1999. 

[14] Christopher M. Kuster, Pierre A. Gremaud, and Rachid Touzani. Fast 
numerical methods for Bernoulli free boundary problems. SIAM J. Sci. 
Comput, 29(2):622-634, 2007. 

[15] R. H. Nochetto, M. Paolini, and C. Verdi. An adaptive finite element 
method for two-phase Stefan problems in two space dimensions. I. Stability 
and error estimates. Math. Comp., 57(195):73-108, Sl-Sll, 1991. 

[16] R. H. Nochetto, M. Paolini, and C. Verdi. An adaptive finite element 
method for two-phase Stefan problems in two space dimensions. II. Im- 
plementation and numerical experiments. SIAM I. Sci. Statist. Comput., 
12(5):1207-1244, 1991. 



17 



[17] Eduardo V. Teixeira Raimundo Leito, Olivaine S. de Queiroz. Regularity 
for degenerate two-phase free boundary problems, eprint arXiv:1202.5264, 
2012. 



[18] L. I. Rubenstein. The Stefan problem. American Mathematical Society, 
Providence, R.I., 1971. Translated from the Russian by A. D. Solomon, 
Translations of Mathematical Monographs, Vol. 27. 

[19] Patricia Saavedra and L. Ridgway Scott. Variational formulation of a model 
free-boundary problem. Math. Comp., 57(196):451-475, 1991. 

[20] K. G. van der Zee, E. H. van Brummelen, and R. de Borst. Goal-oriented 
error estimation and adaptivity for free-boundary problems: the domain- 
map linearization approach. SIAMJ. Sci. Comput, 32(2):1064-1092, 2010. 

[21] Augusto Visintin. Introduction to Stefan-type problems. In Handbook of 
differential equations: evolutionary equations. Vol. IV, Handb. Differ. Equ., 
pages 377-484. Elsevier/North- Holland, Amsterdam, 2008. 

[22] Zhimin Zhang and Ivo Babuska. A numerical method for steady state free 
boundary problems. SI AM J. Numer. Anal, 33(6):2184-2214, 1996. 



18 



